
% HHVR_EM1a_3_types.m is linked to HHVR_EM1b_3_types. 
% In  this example of EM algorithm
% we have three kinds of voters: 
% expressive voters (_u)
% expected utility maximizers (_e)
% frontrunner voters (_f)


function  [parameters_EM   z Z  LogLik loglik   PI total_vote_type ] = HHVR_EM1a_3_types(initial_guess ) % sometimes include ze_fe mean_ze_fe me ze

global DATA choices y y_short X M zz individuals valid_rounds
iterations=333;
convergence =  0.4995;
s=size(DATA,1)/choices;
if s==2304 % i.e. M==1, it behaves badly at the very end, very close to convergence (unknown reason). Does NOT change results.  
    convergence = 5;
end


 
LogLik = [2222 1 ] ;
bbeta = initial_guess ; 
pi_u = 1/3; pi_e = 1/3;pi_f = 1/3; 
s=size(DATA,1)/choices;     % individuals x choices
nume1=zeros(s,1);nume2=zeros(s,1);nume3=zeros(s,1) ;nume4=zeros(s,1) ;nume5=zeros(s,1) ; 
numu1=zeros(s,1);numu2=zeros(s,1);numu3=zeros(s,1) ;numu4=zeros(s,1) ;numu5=zeros(s,1) ;
numf1=zeros(s,1);numf2=zeros(s,1);numf3=zeros(s,1) ;numf4=zeros(s,1) ;numf5=zeros(s,1) ;
prob_e1=zeros(s,1);prob_e2=zeros(s,1);prob_e3=zeros(s,1); prob_e4=zeros(s,1);prob_e5=zeros(s,1); 
prob_u1=zeros(s,1);prob_u2=zeros(s,1);prob_u3=zeros(s,1); prob_u4=zeros(s,1);prob_u5=zeros(s,1); 
prob_f1=zeros(s,1);prob_f2=zeros(s,1);prob_f3=zeros(s,1); prob_f4=zeros(s,1);prob_f5=zeros(s,1); 

denu= zeros(s,1); dene= zeros(s,1);   denf= zeros(s,1); 
denu_x1= zeros(s,1); dene_x1= zeros(s,1); denf_x1= zeros(s,1);


Z=[];iter=0;



s=size(DATA,1)/choices;     % individuals x choices
% numc1=zeros(s,1);numc2=zeros(s,1);numc3=zeros(s,1) ;numc4=zeros(s,1) ;numc5=zeros(s,1) ;numc6=zeros(s,1) ;
% numc7=zeros(s,1) ; 
% nump1=zeros(s,1);nump2=zeros(s,1);nump3=zeros(s,1); nump4=zeros(s,1);nump5=zeros(s,1);nump6=zeros(s,1); 
% nump7=zeros(s,1); prob_c1=zeros(s,1);prob_c2=zeros(s,1);prob_c3=zeros(s,1); prob_c4=zeros(s,1);prob_c5=zeros(s,1);prob_c6=zeros(s,1); 
% prob_c7=zeros(s,1);   
% prob_p1=zeros(s,1);prob_p2=zeros(s,1);prob_p3=zeros(s,1); prob_p4=zeros(s,1); 
% prob_p5=zeros(s,1); prob_p6=zeros(s,1); prob_p7=zeros(s,1);prob_p8=zeros(s,1);  
% denc= zeros(s,1); denp= zeros(s,1); den_cx1= zeros(s,1);denp_x1= zeros(s,1);

Z=[];
for t=1:iterations
            while abs(LogLik(end-1)-LogLik(end)) > convergence
                   DIFF_LOG_LIK= (LogLik(end-1)-LogLik(end) ) 
for i=1:(size(DATA,1))/choices                    % ie from 1 to number of individuals
       
     nume1 (i) =  exp(X(i*choices-4,1)*bbeta(5)    ) ;
     nume2 (i) =  exp(X(i*choices-3,1)*bbeta(5)  + 1*bbeta(1)  ) ;
     nume3 (i) =  exp(X(i*choices-2,1)*bbeta(5)  + 1*bbeta(2)  ) ;
     nume4 (i) =  exp(X(i*choices-1,1)*bbeta(5)  + 1*bbeta(3)  ) ;
     nume5 (i) =  exp(X(i*choices  ,1)*bbeta(5)  + 1*bbeta(4)  ) ;
     dene(i)   =  nume1 (i)  +  nume2 (i) + nume3 (i)+ nume4 (i)+ nume5 (i) ;
     dene_x1(i)   =  nume1(i)*X(i*choices-4,1)  +  nume2(i)*X(i*choices-3,1) + nume3(i)*X(i*choices-2,1)+ ... 
                     nume4(i)*X(i*choices-1,1) +   nume5(i)*X(i*choices  ,1 ) ;
     prob_e1(i)     =  nume1 (i) /  dene(i)  ;     prob_e2(i)     =  nume2 (i) /  dene(i)  ;
     prob_e3(i)     =  nume3 (i) /  dene(i)  ;     prob_e4(i)     =  nume4 (i) /  dene(i)  ;
     prob_e5(i)     =  nume5 (i) /  dene(i)  ;  
    
     numu1 (i) =  exp(X(i*choices-4,2)*bbeta(6)    ) ;
     numu2 (i) =  exp(X(i*choices-3,2)*bbeta(6)  + 1*bbeta(1)  ) ;
     numu3 (i) =  exp(X(i*choices-2,2)*bbeta(6)  + 1*bbeta(2)  ) ;
     numu4 (i) =  exp(X(i*choices-1,2)*bbeta(6)  + 1*bbeta(3)  ) ;
     numu5 (i) =  exp(X(i*choices  ,2)*bbeta(6)  + 1*bbeta(4)  ) ;
     denu(i)   =  numu1 (i)  +  numu2 (i) + numu3 (i)+ numu4 (i)+ numu5 (i) ;
     denu_x1(i)   =  numu1(i)*X(i*choices-4,2)  +  numu2(i)*X(i*choices-3,2) + numu3(i)*X(i*choices-2,2)+ ... 
                     numu4(i)*X(i*choices-1,2) +   numu5(i)*X(i*choices  ,2 ) ;
     prob_u1(i)     =  numu1 (i) /  denu(i)  ;     prob_u2(i)     =  numu2 (i) /  denu(i)  ;
     prob_u3(i)     =  numu3 (i) /  denu(i)  ;     prob_u4(i)     =  numu4 (i) /  denu(i)  ;
     prob_u5(i)     =  numu5 (i) /  denu(i)  ;  

     numf1 (i) =  exp(X(i*choices-4,3)*bbeta(7)    ) ;
     numf2 (i) =  exp(X(i*choices-3,3)*bbeta(7)  + 1*bbeta(1)  ) ;
     numf3 (i) =  exp(X(i*choices-2,3)*bbeta(7)  + 1*bbeta(2)  ) ;
     numf4 (i) =  exp(X(i*choices-1,3)*bbeta(7)  + 1*bbeta(3)  ) ;
     numf5 (i) =  exp(X(i*choices  ,3)*bbeta(7)  + 1*bbeta(4)  ) ;
     denf(i)   =  numf1 (i)  +  numf2 (i) + numf3 (i)+ numf4 (i)+ numf5 (i) ;
     denf_x1(i)   =  numf1(i)*X(i*choices-4,3) +   numf2(i)*X(i*choices-3,3) + numf3(i)*X(i*choices-2,3)+ ... 
                     numf4(i)*X(i*choices-1,3) +   numf5(i)*X(i*choices  ,3 ) ;
     prob_f1(i)     =  numf1 (i) /  denf(i)  ;     prob_f2(i)     =  numf2 (i) /  denf(i)  ;
     prob_f3(i)     =  numf3 (i) /  denf(i)  ;     prob_f4(i)     =  numf4 (i) /  denf(i)  ;
     prob_f5(i)     =  numf5 (i) /  denf(i)  ; 

end
 
densitye=zeros(s,1)'; densityu=zeros(s,1)'; densityf=zeros(s,1)';
        
 
for i=1:(size(DATA,1))/choices
    if y_short(i)==1
        densitye(i)= prob_e1(i)  ;
        densityu(i)= prob_u1(i) ;
        densityf(i)= prob_f1(i) ;
    elseif y_short(i)==2
        densitye(i)= prob_e2(i) ;
        densityu(i)= prob_u2(i) ;
        densityf(i)= prob_f2(i) ;
    elseif y_short(i)==3
        densitye(i)= prob_e3(i) ;
        densityu(i)= prob_u3(i) ;
        densityf(i)= prob_f3(i) ;
    elseif y_short(i)==4
        densitye(i)= prob_e4(i) ;
        densityu(i)= prob_u4(i) ;
        densityf(i)= prob_f4(i) ;
    elseif y_short(i)==5
        densitye(i)= prob_e5(i) ;
        densityu(i)= prob_u5(i) ;
        densityf(i)= prob_f5(i) ;
     end
end

density=[densitye ; densityu ; densityf] ;         % 3  x individuals and they are all different
ze=zeros((size(DATA,1))/choices,1);                % individuals  x  1
zu=zeros((size(DATA,1))/choices,1);                % individuals  x  1
zf=zeros((size(DATA,1))/choices,1);                % individuals  x  1

for i=1:(size(DATA,1))/choices
    ze(i) = ( pi_e*densitye(i)) / ( pi_e*densitye(i) + pi_u*densityu(i) +pi_f*densityf(i)  );
    zu(i) = ( pi_u*densityu(i)) / ( pi_e*densitye(i) + pi_u*densityu(i) +pi_f*densityf(i)  );
    zf(i) = ( pi_f*densityf(i)) / ( pi_e*densitye(i) + pi_u*densityu(i) +pi_f*densityf(i)  );
 end
Z=[Z ze zu zf];
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % FE: averaging the probabilities for each individual
%  
%     ze_fe = reshape(ze,  valid_rounds, individuals/valid_rounds);   % valid_rounds(48) x individuals   
%     zu_fe = reshape(zu,  valid_rounds, individuals/valid_rounds);
%     zf_fe = reshape(zf,  valid_rounds, individuals/valid_rounds);
%  
%   % Next: means for each individual (taken from individual-action)
% mean_ze_fe = mean(ze_fe); mean_zu_fe = mean(zu_fe);mean_zf_fe = mean(zf_fe);
% me = repmat(mean_ze_fe,   ( valid_rounds),1 );                     % valid_rounds(48) x individuals 
% mu = repmat(mean_zu_fe,   ( valid_rounds),1 );
% mf = repmat(mean_zf_fe,   ( valid_rounds),1 );
% ze = reshape(me, numel(me),1);zu = reshape(mu, numel(me),1); zf = reshape(mf, numel(me),1);
   
  % FE creation ends up  here 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
 
            loglik=[];Loglikvuong = [];
            for i=1:s
                loglik(i) =  ze(i)*log(densitye(i)) + zu(i)*log(densityu(i)) + zf(i)*log(densityf(i))  ;
                Loglikvuong = [Loglikvuong ; loglik(i)];
            end
            Loglik = sum(loglik);
            LogLik = [LogLik Loglik]  ;
            DIFF_LOG_LIK= (LogLik(end-1)-LogLik(end) );
            pi_e = mean(ze); pi_u = mean(zu); pi_f = mean(zf); 
            LogLik(end);
            PI = [pi_e pi_u pi_f ] ;
            
zz = [ze zu zf ] ;
inputs = [bbeta(1:7) ]; 
 
[bbeta ] = fsolve(@(bbeta)    HHVR_EM1b_3_types(bbeta), inputs,optimset('maxiter',20000       , 'MaxFunEvals',40000      ) ) ; % not [bbeta,fval,a,b]
 
CONSTANT =  [bbeta(1) bbeta(2)  bbeta(3) bbeta(4)    ];
BBETA = [bbeta(5) bbeta(6)  bbeta(7)     ]; 
iter=iter+1;

            end %while
end

if size(DATA,2)>12
id_agents = reshape(DATA(:,1)',5,numel(DATA(:,1))/5) ;
id_agents = id_agents(1,:);
disjoint = reshape(DATA(:,13)',5,numel(DATA(:,13))/5) ;
disjoint = disjoint(1,:);
size(disjoint)
LogLik =  LogLik';
z=[ zz id_agents' disjoint'];
else
    z = zz;
end
 
total_vote_type = [0 0 0 ];
for i = 1: size(DATA,1)
    if DATA(i,2)== 1 && DATA(i,3) ==1
        total_vote_type(1) = total_vote_type(1) + 1;
    elseif DATA(i,2)== 1 && DATA(i,4) ==1
        total_vote_type(2) = total_vote_type(2) + 1;
    elseif DATA(i,2)== 1 && DATA(i,5) ==1
        total_vote_type(3) = total_vote_type(3) + 1;
    end
end
parameters_EM = [bbeta ];



 



